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MULTI-DIMENSIONAL SUMMATION-BY-PARTS OPERATORS: 
GENERAL THEORY AND APPLICATION TO SIMPLEX ELEMENTS 
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Abstract. Summation-by-parts (SBP) finite-difference discretizations share many attractive 
properties with Galerkin finite-element methods (FEMs), including time stability and superconver- 
gent functionals; however, unlike FEMs, SBP operators are not completely determined by a basis, so 
the potential exists to tailor SBP operators to meet different objectives. To date, application of high- 
order SBP discretizations to multiple dimensions has been limited to tensor product domains. This 
paper presents a definition for multi-dimensional SBP finite-difference operators that is a natural ex¬ 
tension of one-dimensional SBP operators. Theoretical implications of the definition are investigated 
for the special case of a diagonal norm (mass) matrix. In particular, a diagonal-norm SBP operator 
exists on a given domain if and only if there is a cubature rule with positive weights on that domain 
and the polynomial-basis matrix has full rank when evaluated at the cubature nodes. Appropriate 
simultaneous-approximation terms are developed to impose boundary conditions weakly, and the 
resulting discretizations are shown to be time stable. Concrete examples of multi-dimensional SBP 
operators are constructed for the triangle and tetrahedron; similarities and differences with spectral- 
element and spectral-difference methods are discussed. An assembly process is described that builds 
diagonal-norm SBP operators on a global domain from element-level operators. Numerical results 
of linear advection on a doubly periodic domain demonstrate the accuracy and time stability of the 
simplex operators. 

Key words, summation-by-parts, finite-difference method, unstructured grid, spectral-element 
method, spectral-difference method, mimetic discretization 
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1. Introduction. Summation-by-parts (SBP) operators are high-order finite- 
difference schemes that mimic the symmetry properties of the differential operators 
they approximate [19]. Respecting such symmetries has important implications; in 
particular, they enable SBP discretizations that are both time stable and high-order 
accurate [4, 34, 27], properties that are essential for robust, long-time simulations of 
turbulent flows [25, 35]. 

Most existing SBP operators are one-dimensional [30, 24, 31, 23] and are applied 
to multi-dimensional problems using a multi-block tensor-product formulation [32, 14, 
29]. Like other tensor-product methods, the restriction to multi-block grids compli¬ 
cates mesh generation and adaptation, and it limits the geometric complexity that 
can be considered in practice. 

The limitations of the tensor-product formulation motivate our interest in gener¬ 
alizing SBP operators to unstructured grids. There are two ways this generalization 
has been pursued in the literature: 1) construct global SBP operators on an arbi¬ 
trary distribution of nodes, or; 2) construct SBP operators on reference elements and 
assemble a global discretization by coupling these smaller elements. 

The first approach is appealing conceptually, and it is certainly viable for second- 
order accurate SBP schemes. For example, Nordstrom et al. [28] showed that the 
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vertex-centered second-order-accurate finite-volume scheme^ has a multi-dimensional 
SBP property, even on unstructured grids; however, the first approach presents chal¬ 
lenges when constructing high-order operators. Kitson et al. [17] showed that, for 
a given stencil width and design accuracy, there exist grids for which no stable, 
diagonal-norm SBP operator exists. Thus, building stable high-order SBP opera¬ 
tors on arbitrary node distributions may require unacceptably large stencils. When 
SBP operators do exist for a given node distribution, they must be determined glob¬ 
ally by solving a system of equations, in general. The global nature of these SBP 
operators is exemplified in the mesh-free framework of Chiu et al. [20, 5]. 

The second approach — constructing SBP operators on reference elements and 
using these to build the global discretization — is more common and presents fewer dif¬ 
ficulties. The primary challenge here is to extend the one-dimensional SBP operators 
of Kreiss and Scherer [19] to a broader set of operators and domains. The existence 
of such operators, at least in the dense-norm case^, was established by Carpenter and 
Gottlieb [2]. They proved that operators with the SBP property can be constructed 
from the Lagrangian interpolant on nearly arbitrary nodal distributions, which is 
practically feasible on reference elements with relatively few nodes. More recently, 
Gassner [10] showed that the discontinuous spectral-element method is equivalent to 
a diagonal-norm SBP discretization when the Legendre-Gauss-Lobatto nodes are used 
with a lumped mass matrix. 

Of particular relevance to the present work is the extension of the SBP concept 
by Del Rey Fernandez et al. [7]. They introduced a generalized summation-by-parts 
(GSBP) definition for arbitrary node distributions on one-dimensional elements, and 
these ideas helped shape the definition of SBP operators presented herein. 

Our first objective in the present work is to develop a suitable definition for 
multi-dimensional SBP operators on arbitrary grids and to characterize the resulting 
operators theoretically. We note that the discrete-derivative operator presented in [20] 
is a possible candidate for defining (diagonal-norm) multi-dimensional SBP operators; 
however, it lacks properties of conventional SBP operators that we would like to retain, 
such as the accuracy of the discrete divergence theorem [15]. 

Our second objective is to provide a concrete example of multi-dimensional 
diagonal-norm SBP operators on non-tensor-product domains. We focus on diagonal- 
norm operators, because they are better suited to discretizations that conserve non¬ 
quadratic invariants [17]; they are also more attractive than dense norms for explicit 
time-marching schemes. We construct diagonal-norm SBP operators for triangular 
and tetrahedral elements. The resulting operators are similar to those used in the 
nodal triangular-spectral-element method [6, 26, 11]. Unlike the spectral-element 
method based on cubature points, the SBP method is not completely specified by 
a polynomial basis; we use the resulting freedom to enforce the summation-by-parts 
property, which leads to provably time-stable schemes. 

The remaining paper is structured as follows. Section 2 presents notation and the 
proposed definition for multi-dimensional SBP operators. We study the theoretical 
implications of the proposed definition in Section 3. We then describe, in Section 4, 
how to construct diagonal-norm SBP operators for the triangle and tetrahedron. Sec¬ 
tion 4 also establishes that SBP operators on subdomains can be assembled into SBP 
operators on the global domain. Results of applying the triangular SBP operators 


^On simplices, the vertex-centered finite-volume scheme is equivalent to a mass-lumped p = 1 
finite-element discretization 
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to the linear advection equation are presented in Section 5. Conclusions are given in 
Section 6. 

2. Preliminaries. To make the presentation concise, we concentrate on multi¬ 
dimensional SBP operators in two dimensions; the extension to higher dimensions 
follows in a straightforward manner. Furthermore, we present definitions and theo¬ 
rems for operators in the x coordinate direction only; the corresponding definitions 
and theorems for the y coordinate direction follow directly from those in the x direc¬ 
tion. 

2.1. Notation. We consider discretized derivative operators defined on a set of n 
nodes, S = {(a^i, 2/i)}r=i- Capital letters with script type are used to denote continuous 
functions. For example, U{x) € L^(0) denotes a square-integrable function on the 
domain 11. We use lower-case bold font to denote the restriction of functions to the 
nodes. Thus, the restriction of W to S' is given by 

u = \lA{xi,yi ),. . . ,U{Xn,yn)\^ ■ 

Several theorems and proofs make use of the monomial basis. For two spatial 
variables, the size of the polynomial basis of total degree p is 

_ (P+ l)(P + 2) 

P- 2 

More generally, n* = where d is the spatial dimension. We use the following 

single-subscript notation for monomial basis functions: 

Vk{x,y) = x'‘y^~\ k= j{j + l)/2 + i + l, V j e {0,1,... ,p}, z e {0,1,..., j}. 

We will frequently evaluate Vk and dVk/dx at the nodes S, so we introduce the 
notation 


Pk = [Vk{xi,yi), . . . ,Vk{Xn,yn)V > 

ap,, , crPt, 

Finally, matrices are represented using capital letters with sans-serif font; for 
example, the first-derivative operators with respect to x and y are represented by the 
matrices and Dy, respectively. Entries of a matrix are indicated with subscripts, 
and we follow Matlab®-like notation when referencing submatrices. For example, A. j 
denotes the column of matrix A, and A. i-j, denotes its first k columns. 

2.2. Multi-dimensional SBP operator definition. We propose the following 
definition for D^,, the SBP first-derivative operator with respect to x. An analogous 
definition holds for and, in three-dimensions, D^. Definition 2.1 is a natural 
extension of the definition of GSBP operators proposed in [7], which itself extends 
the classical SBP operators introduced by Kreiss and Scherer [19]. 

Definition 2.1. Two-dimensional summation-by-parts operator: Con¬ 
sider an open and bounded domain D C with a piecewise-smooth boundary F. The 
matrix Dx is a degree p SBP approximation to the first derivative on the nodes 

S = {(a:i,2/i)}r=i */ 

1. DxPk=p'k, Vfc e {l,2,...,n*},- 

2. Dj, = where H is symmetric positive-definite, and; 


and p'j, = 
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3. Sj; = Qj, + where Qj = —Q^, = E^,, and satisfies 

Pk ^xPm V /c, TTl G {1,2,..., , 


i 


where t > p and Ux is the x component of n = [ux^ny^ , the outward pointing unit 


normal on F. 

Before studying the implications of Definition 2.1 in Section 3, it is worthwhile 
to motivate and elaborate on the three properties in the definition. 

Property 1 ensures that D^, is an accurate approximation to the first partial 
derivative with respect to x. The operator must be exact for polynomials of total 
degree less than or equal to p, so at least n* nodes are necessary to satisfy property 1 . 

Remark 1. We emphasize that a polynomial basis is not used to define the 
solution in SBP methods, in contrast with the piecewise polynomial expansions found 
in finite-element methods. We adopt the monomial basis only to define the accuracy 
conditions concisely and avoid cumbersome Taylor-series expansions. 

The matrix H must be symmetric positive-definite to guarantee stability: without 
property 2, the discrete “energy”, u^Hu, could be negative when vf^u > 0, and vice 
versa. The so-called norm matrix H can be interpreted as a mass matrix, i.e. 



but it is important to emphasize that SBP operators are finite-difference operators, 
and there is no (known) closed-form expression for an SBP nodal basis in 

general. In the diagonal norm case, we shall show that another interpretation of H is 
as a cubature rule. 

Property 3 is needed to mimic integration by parts (IBP). Recall that the IBP 
formula for the x derivative is 



The discrete version of the IBP formula, which follows from property 3, is 
v^UDxU-\-vfi^UDxV = v^ExU, 

There is a one-to-one correspondence between each term in the IBP formula and its 
SBP proxy. For example, it is clear from property 3 that v^ExU approximates the 
surface integral in IBP to order r. Moreover, in Section 3 we show that diagonal-norm 
SBP operators also approximate the left-hand side of IBP. 

3. Analysis of diagonal-norm, multi-dimensional, summation-by-parts 
operators. In this section, we determine the implications of Definition 2.1 on the 
constituent matrices of a multi-dimensional SBP operator and whether or not such 
operators exist. We also investigate the time stability of discretizations based on 
multi-dimensional SBP operators. The focus is on diagonal-norm operators; however, 
the ideas presented here can be extended to dense-norm operators, i.e. where the 
matrix H is not diagonal. 

The following lemma will prove useful in the sequel. It follows immediately from 
properties 1 and 3, so we state it without proof. 

Lemma 3.1 (compatibility). Let D^, = H~^(Qa; -|- ^E^,) be an SBP operator of 
degree p. Then we have the following set of relations: 


pl,^p'k+Pk^P'm^Pli^xPk, \/k,mG {l,2,...,n*p}. 


(3.1) 
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We refer to (3.1) as the compatibility equations for the x derivative; H must 
simultaneously satisfy analogous relations for Ey. The relation between H and E^; 
was first derived by Kreiss and Scherer [19] and Strand [30] to construct a theory 
for one-dimensional classical hnite-difference-SBP operators. Furthermore, Del Rey 
Fernandez et al. [7] have used these relations to extend the theory of such operators to 
a broader set. What is presented in this paper is a natural extension of those works to 
multi-dimensional operators, and the derivation of (3.1) follows in a straightforward 
manner from any of the mentioned works. 

Our first use of the compatibility equations is to prove that, in the diagonal-norm 
case, the multi-dimensional SBP definition conceals a cubature rule with positive 
weights. 

Theorem 3.2. Let = H~^Sa; he a degree p, diagonal-norm, multi-dimensional 
SBP operator on the domain Q. Then the nodes S = {(a^i,?/i)}?=i diagonal 

entries of H form a degree 2p — 1 cubature rule on Q. 

Proof. Using property 3 of Definition 2.1, the compatibility equations become 

n 

EAj 




dx 


dx 


= 4> VmVknxdV, 


V fc, TO G {1,2,..., n* 


Using the chain rule on the left and integration by parts on the right results in 




d'Pm'Pk 

dx 




/ 


d'Pm'Pk 

dx 


dO, 


yk,me { 1 , 2 ,..., 71 *}. 


(3.2) 


Since Pk and Pm are monomials of degree at most p, it follows that d {PmPi) /dx is a 
scaled monomial of degree at most 2p— 1; thus, by considering all of the combinations 
of k and to, (3.2) implies 


HjjPfe {xj, Uj) 

1=1 


Jpkdn, vfcG{i,2,...,77;p_j, 
a 


which are the conditions for a degree 2p — 1 cubature. □ 

We now prove one of our central theoretical results, relating the existence of a 
diagonal-norm SBP operator to the existence of a cubature rule with positive weights. 

Theorem 3.3. Consider the node set S = {{xi,yi)}'^^i with n > n* nodes, 
and define the generalized Vandermonde matrix V G whose columns are the 

monomial basis evaluated at the nodes; 


y:,k=Pk, Vfc G {1,2, ...,n*} 

Assume that the columns of V are linearly independent. Then the existence of a 
cubature rule of degree at least 2p — 1 with positive weights on S is necessary and 
sufficient for the existence of degree p diagonal-norm SBP operators approximating 
the first derivatives and on the node set S. 

Proof. The necessary condition on H follows immediately from Theorem 3.2. To 
prove sufficiency, we must show that, given a cubature rule, we can construct an 
operator that satisfies properties 1-3 of Definition 2.1 on the same node set as the 
cubature rule. 
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Before proceeding, we introduce some matrices that facilitate the proof. Let 
Vj, G be the matrix whose columns are the x derivatives of the monomial- 

basis: 

(V-),k=Pk, Vk€{l,2,...,n;}. 

We construct an invertible matrix V € by appending a set of vectors, W G 

jgnx(n-np) ^ that are linearly independent amongst themselves and to the vectors in 
V: 


V = [V w]. 

Similarly, we define 

V, = [V, w,], 

where G is matrix that will be specified later. Below, we use the 

degrees of freedom in to satisfy the SBP definition. 

Let H be the diagonal matrix whose entries are the cubature weights ordered 
consistently with the cubature node set S. Since the cubature weights are positive, 
property 2 is satisfied. 

Next, we use the cubature to construct a suitable Ea,. Using V and Vx, we define 
the symmetric matrix 

Ea = V^HV. + VjHV. (3.3) 


Since V and are polynomials of degree p and p — 1, respectively, evaluated at the 
nodes, the cubature is exact for the right-hand side of (3.3): 


(^■), 


+ / n^cin = / PiP^n^dr, 

OX Jo ox Jr 


(3.4) 


V fc, m £ {1, 2,..., n*}. Now we can define the boundary operator 


- w-T 


E. = V 


E. FT 

F. t G.T 


V' 


where G and = G]! G r("-’^;)x(’^-";). It follows from this definition 

that Ex is symmetric. Moreover, together with (3.4), this definition implies 

(VTe,V)^^ = (e,) = (f Vk'PmnxdT, V A:,m G {1, 2,..., n;}, 

’ \ / fc,m Jy' 


so Ea; satisfies the accuracy condition of property 3. 
Finally, let 

Qx = HV,V-i - ^E,. 


(3.5) 


The accuracy conditions, which are equivalent to showing Da^V = Vx, follow immedi¬ 
ately from this definition of Qa;: 

DaV = + ^E,^ V = H-i (hV,V-i) V = V,; 
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thus, property 1 is satisfied. 

Our remaining task is to show that Q^, can be constructed to be antisymmetric. 
If we can show that 


VTQ.V VTQ.VV' 
wtq^v \N'^Q^\N 

can be made antisymmetric, then the result will follow for Q^. Consider the first 
block in the 2x2 block matrix above, i.e. 

V'^Q,,V = - ^V^E^V. 

Adding this block to its transpose, we find 

V'^Q,,V + V^QjV = V^HV,, -P VjHV - V^E^V, (3.6) 

where we have used the symmetry of Ex- The right-hand side of (3.6) is the matrix 
form of the (rearranged) compatibility equations (3.1). Thus, V'^Q 3 ,V + = 

0, proving that the first block is antisymmetric. For the remaining three blocks, 
antisymmetry requires 

(V^Q^W)^ =-W^Q^V, and =-W'^QJW. 

Substituting Qx and Ex and simplifying, we obtain the following equations: 

WjHV-1-W'^HV^ = and W^HW^-F WJHW = G^. (3.7) 

The two matrix equations above constitute n(n — n*) scalar equations. We are free 
to choose Wj,, Gxi and F^,, so the matrix equations are underdetermined (W^, alone 
has n{n — n*) entries). To prove existence of an SBP operator we need only find one 
solution; for example, take W^, = 0, G^, = 0, and Ex = W^^HVa,. Thus, the equations 
can be satisfied to ensure the antisymmetry of Qa,. □ 

Remark 2. In general, there are infinitely many operators associated with a 
given cubature rule that satisfy Definition 2.1. For example, the proof of Theorem 3.3 
only considered one way to solve (3.7). Another way to satisfy these conditions is to 
set Fa: = 0 and then solve 


V^QxV = 


(V^H) Wa: = -VJHW 

for \Nx by finding the minimum Frobenius-norm solution. Gx can then be computed 
directly from the second eguation, (3.7). 

The following theorem characterizes the matrices Sa; and Qa; in terms of the bilin¬ 
ear forms that they approximate. The theorem is useful when discretizing the weak 
form of a PDF, rather than the strong form, using SBP finite-difference operators. 

Theorem 3.4. Let Da; = H“^Sa: = (Qa; + be a diagonal-norm SBP 
operator of degree p on the domain fl. Then 

Pk^xPm = j Vk^^dD, V fe-F m < min (n 2 p,n*) (3.8) 

plQxPm= [ Vk^^dD-^ 1 VkVmnxdT, V fc -F m < min (n^p, n*) . (3.9) 
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Proof. The proof of (3.8) is analogous to the proof in [15] for one-dimensional 
classical hnite-difference-SBP operators. With (3.8) established, we can substitute 
Sx = Qa; + 5 Ea; and rearrange to obtain 


Pfc QxPm 


f T, 1 Tp 


V fc -|- TO < min (n 2 p, n*) , 


where we have used the accuracy property of E^; to get the desired result. 0 


3.1. Stability Analysis. We conclude our analysis of diagonal-norm multi¬ 
dimensional SBP operators by investigating the stability of an SBP semi-discretization 
of the constant-coefficient advection equation 


dt 


' dx 


3^=0 

^dy 


U{x,y,t) = Uhc{x,y,t), 


V {x,y) € O, 
V {x,y) G P-, 


(3.10) 


where P” = {(a;, y) G P | Pxnx + PyTiy < 0} is the inflow boundary and P"*" = P \ P“ 
is the outflow boundary. 

Let Da; and Dj, be SBP operators and let and Ey be their corresponding 
boundary operators. In order to impose the boundary conditions in a stable manner, 
we introduce the decomposition 


l3xEx + l3yEy — E_|_ -|- E_, (3.11) 

where E+ is symmetric positive semi-definite and E_ is symmetric negative semi- 
definite, and 


pIE± pm = f VkVm {Pxiix + PyTiy) dP, W k,mG {1,2,..., nl}. 


r± 


The existence of the decomposition (3.11) is guaranteed provided Fa, = Fj, = 0; see 
Appendix A for the proof. 


Using Da 
given by 


and Dy and the matrix E_, a consistent semi-discretization of (3.10) is 


du 

dt 


+ l3xE)xU P jdyDyU = aY\ ^E_(M-itbc). 


(3.12) 


The three terms on the left-hand side of (3.12) correspond to the three terms in 
the strong-form of the PDE. The terms on the right-hand side of (3.12) are penalties 
that enforce the boundary conditions weakly using simultaneous-approximation terms 
(SATs) [9, 3]. The boundary data is supplied by the vector Ubc, which must produce 
a sufficiently accurate reconstruction of Uhc along the boundary. Evaluating Uhc 
is adequate when the nonzeros of Ex correspond to nodes that lie on P, such as 
the simplex operators presented below. More generally, a preprocessing step can be 
performed to find a suitable Mbc- 

We now show that (3.12) is time stable. 

Theorem 3.5. Letu be the solution to (3.12) with homogeneous boundary eondi- 
tions and bounded initial condition. Then the norm ||it||H = V u"^Hu is non-inereasing 

ifa> i. 
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Proof. Multiplying the semi-discretization (with itbc = 0) from the left by 
we find 


rp ^ rp rp 

u H— -f pxU SxU + /3yU SyU — au E-U 


dt 


= —+ (2(7 — 


To obtain the last line, we used the definitions of S^, and Sy, as well as the decom¬ 
position (3.11). The first quadratic form on the right is non-positive by definition of 
E_|-. The result follows if cr > i, since E_ is negative semi-definite. □ 

4. Constructing the operators. This section describes how we construct 
diagonal-norm SBP operators for triangles and tetrahedra. The algorithms described 
below have been implemented in the Julia package SummationByParts^. 

4.1. The node coordinates and the norm matrix. Theorem 3.2 tells us 
that the diagonal entries in H are positive weights from a cubature that is exact for 
polynomials of total degree 2p — 1. Thus, our first task is to find cubature rules with 
positive weights for the triangle and tetrahedron. Additionally, we seek rules that use 
as few nodes as possible for a given order of accuracy while respecting the symmetries 
of the triangle and tetrahedron; the former condition is for efficiency while the latter 
condition is to reduce directional biases. 

For the operators considered in this work, we require that cubature nodes 

lie on each boundary facet, where d is the spatial dimension. This requirement on 
the nodes is motivated by the particular form of the E^,, Ey, and Ej, operators that 
we consider below; however. Definition 2.1 does not require a prescribed number of 
boundary nodes, and SBP operators for the 2- and 3-simplex may exist that do not 
have any boundary nodes at all. 

Cubature rules that meet our requirements for triangular elements are presented 
in references [21, 6, 26, 11] in the context of the spectral-element and spectral- 
difference methods. Table 1 summarizes the rules that are adopted for triangular- 
element SBP operators of degree p = 1, 2, 3, and 4. For reference, the node locations 
for the triangular cubature rules are shown in Figure 1. 

To find cubature rules for the tetrahedron, we follow the ideas presented in [11, 
36, 33]. Our procedure is briefly outlined below for completeness, but we make no 
claims regarding the novelty of the cubature rules or our method of finding them. 

We assume that each node belongs to a (possibly degenerate) symmetry orbit [36] . 
As indicated above, we assume that the cubature node set includes p -|- 1 nodes along 
each edge and (p+l)(p-|-2)/2 nodes on each triangular face. For the interior nodes, we 
activate the minimum number of symmetry orbits necessary to satisfy the accuracy 
conditions; these orbits have been identified through trial and error. 

Each symmetry orbit has a cubature weight associated with it, and orbits that 
are non-degenerate are parameterized using one or more barycentric parameters. To¬ 
gether, the orbit parameters and the weights are the degrees of freedom that must be 
determined. They are found by solving the nonlinear accuracy conditions using the 
Levenberg-Marquardt algorithm. The accuracy conditions are implemented using the 
integrals of orthogonal polynomials on the tetrahedron [18, 8]. 

Table 2 summarizes the node sets used for the tetrahedron cubature rules, and 
Figure 2 illustrates the node coordinates. 


®https://github.com/OptimalDesignLab/SummationByParts.jl 
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Table 1 

Active orbits and their node counts for triangular-element operators. The notation Perm indi¬ 
cates that every permutation of the barycentric coordinates is to be considered. Free-node counts are 
decomposed into the product of the number of nodes in the orbit times the number of orbits of that 

type- 


operator degree, p 



orbit name 

barycentric form 

1 

2 

3 

4 

fixed nodes 

vertices 

Perm(l, 0,0) 

3 

3 

3 

3 


mid-edge 

Perm 

— 

3 

— 

3 


centroid 

(1 1 i) 

— 

1 

— 

— 

free nodes 

edge 

Perm (a, 1 — a, 0) 

— 

— 

6x1 

6x1 


S '21 

Perm (a, a, 1 — 2a) 

— 

— 

3 X 1 

3x2 



# free parameters 

— 

— 

2 

3 



# nodes total 

3 

7 

12 

18 






Fig. 1. Node distributions for cubature rules adopted for the SBP operators on triangles. 




p = 3 


p = 4 


Fig. 2. Node distributions for cubature rules adopted for the SBP operators on tetrahedra. 



4.2. The boundary operators. Definition 2.1 implies that the boundary op¬ 
erator Ex satisfies 

= j) UVrix dr 

for all polynomials U and V whose total degree is less than t > p. In particular, if we 
choose U and V to be nodal basis functions on the faces, we can isolate the entries of 
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Table 2 

Active orbits and their node counts for tetrahedral-element operators. See the caption of Table 1 
for an explanation of the notation. 





operator degree, p 


orbit name 

barycentric form 

1 2 

3 

4 

fixed nodes 

vertices 

Perm(l, 0,0, 0) 

4 4 

4 

4 


mid-edge 

Perm (i, i, 0, O) 

— 6 

— 

6 


centroid 

(1 L L 1 ) 

V4’ 4’ 4’ 4/ 

— 1 

— 

1 


face centroid 

Perm(i,i,i,0) 

~ — 

4 

— 

free nodes 

edge 

Perm {a, 1 — 0,0,0) 

— — 

12 X 1 

12 X 1 


face S21 

Perm (a, a,l — 2 a, 0) 

— — 

— 

12 X 1 


S31 

Perm {a, a,a,l — 3 a) 

— — 

4 X 1 

4x1 


S22 

Perm (o, a,\ — a,\ — a) 

— — 

— 

6x1 



ff free parameters 

- - 

2 

4 



nodes total 

4 11 

24 

45 


£x- This is possible because we have insisted on operators with nodes on each 

facet, which leads to a complete nodal basis. For further details on the construction 
of the boundary operators, we direct the interested reader to [13, pg. 187]. 

Remark 3. The boundary operators, when restricted to the boundary nodes, 
are dense matrices. Contrast this with the tensor-product case, where the boundary 
operators are diagonal matrices. In the simplex case, we have not found a way to 
construct diagonal E^, Ey and that are sufficiently accurate. 

4.3. The antisymmetric part. The accuracy conditions are used to determine 
the antisymmetric matrices Qx, Qy, and Q^. We will describe the process for Q^,, since 
it can be adapted in a straightforward way to Qy and, in the case of the tetrahedron, 

Qz. 

In theory, we can compute Qx using the monomials that appear in the SBP 
operator definition; however, the monomials are known to produce ill-conditioned 
Vandermonde matrices. Instead, we follow the standard practice in spectral-element 
methods and apply the accuracy conditions to appropriate orthogonal bases on the 
triangle and tetrahedron [18, 8, 13]. Unlike finite- and spectral-element methods, the 
basis alone does not completely specify an SBP operator. 

Let P and P^, be the matrices whose columns are the orthogonal basis function 
values and derivatives, respectively, evaluated at the nodes. Then the accuracy con¬ 
ditions imply Da;P = P^,, or, in terms of the unknown Qx, 

QxP = HP,-iE,P. 

This can be recast as the linear system 


Aq = b, 


(4.1) 
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where q denotes a vector whose entries are the strictly lower part of Q^,: 

g( 0 - 2 Ki-i) ^ , 2<i<n,l<j <i. 

There are (n — l)n/2 unknowns and n x equations in (4.1); thus, for the op¬ 

erators considered here, there are more equations than unknowns. Fortunately, the 
compatibility conditions ensure that the system is consistent. Indeed, for p > 3 the 
rank of A is actually less than the size of q, so there are an infinite number of solutions. 
In these cases, we choose the minimum-norm least-squares solution [12] to (4.1). 

4.4. Similarities and differences with existing operators. There is a vast 
literature on high-order discretizations for simplex elements, so we focus on the two 
that share the most in common with the proposed SBP operators: the diagonal 
mass-matrix spectral-element (SE) method [6, 26, 11] and the spectral-difference (SD) 
method [22]. 

Our norm matrices H are identical to the lumped mass matrices in the SE method. 
The difference between the methods arises in the definition of S^. In the SE method 
of Giraldo and Taylor [11], the S^, matrix is defined as 

^ dcj) 

i^x ) i j ~ ^ ^ (xfc , J/fc) {^ki yk')i 

k=0 

where {(l>i}^=i is the so-called cardinal basis. For a degree p operator, the cardinal 
basis is a polynomial nodal basis that is a super-set of the basis for degree p polyno¬ 
mials; the basis contains polynomials of degree greater than p, because the number 
of nodes n is greater than , in general. Consequently, the cubature defined by 
H is not exact for the product (j)id(j)j/dx when p > 2, and the resulting S®® does 
not satisfy the SBP definition for the p > 2 discretizations. Indeed, as the results 
below demonstrate, the higher-order SE operators are unstable and require filtering 
or numerical dissipation even for linear problems. 

Remark 4. The matrix in the diagonal mass-matrix SE method can be made 
to satisfy the SBP definition by using a cubature rule that is exaet for the eardinal 
basis; however, such a cubature rule would require additional cubature points and would 
defeat the purpose (i.e. efficiency) of collocating the cubature and basis nodes. 

Diagonal-norm SBP operators can also be viewed as a special case of the SD 
method in which the unknowns and fluxes are collocated. As pointed out in [22] , this 
means that our proposed SBP operators require more unknowns to achieve a given 
accuracy than spectral-difference methods; however, collocation eliminates the recon¬ 
struction step, so there is a tradeoff between memory and floating-point operations. 
More importantly, this relative increase in unknowns applies only to discontinuous 
methods. If we assemble a global SBP operator, as described below, then the number 
of unknowns can be significantly reduced. 

4.5. Assembly of global SBP operators from elemental operators. The 

SBP operators defined in Sections 4.1-4.3 can be used in a nodal DC formula¬ 
tion [13] with elements coupled weakly using, for example, simultaneous approxima¬ 
tion terms [9, 3]. An alternative use for these element-based operators, and the one 
pursued here, is to mimic the continuous Galerkin formulation. That is, we assemble 
global SBP operators from the elemental ones. 

We need to introduce some additional notation to help describe the assembly 
process and facilitate the proof of Theorem 4.1 below. Suppose the domain D is 
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partitioned into a set of L non-overlapping subdomains with boundaries 

L 

= and n O(') = 0, V fc /, 

1=1 

where U denotes the closure of 

Each subdomain is associated with a set of nodes \such 

that G In the present context, some of the nodes in will lie on the 

boundary F*^*) and be shared by adjacent subdomains. 

Let S = Suppose there are n unique nodes in S', and let each node be 

assigned a unique global index. Suppose i is the global index corresponding to the 
local node of element 1. We define 7S^\i,j) to be the n x h matrix with zeros 
everywhere except in the (q j) entry, which is unity. If ej denotes the i column of the 
n X h identity, then = e^ej. 

We can now state and prove the following result. 

Theorem 4.1. Let ^ be a degree p SBP operator for the first 

derivative d/dx on the node set If 


L n 


L n n 


1=1 i=i j=i 


then Da; = H~^Sa; is a degree p SBP operator on the global node set S. 

Proof. We need to check each of the three properties in Definition 2.1. 

1. The first property is straightforward, albeit tedious, to verify. exists by 
property 2, which is shown to hold independently below, so we have 


DxPk = H Sa;Pk 

L n n 

= H-i ' 


1 _1 _1 X_1 


1=1 i=i j=i 

jW 


Pk 


H ^EE ,,(0 E eiT’fe(a;j,yj) 

1^1 2=1 j = l 


L n 




1=1 i=l 
L 


i=i 




1=1 i=l 
L n 


= H- 




1=1 i=l 


Pk, 


But the term in brackets above is the definition of H, so we are left with Da;Pfc = pj-, 
as desired. 

2. H is clearly diagonal and positive by construction, so property 2 is satisfied. 
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3. Finally, we form the decomposition where 


L n n 

L n n 

/=! i=l 3 = 1 


The matrix is antisymmetric, because it is the sum of antisymmetric matrices. 
Similarly, is symmetric, because it is the sum of symmetric matrices. In addition. 


L n n 

pJ^xPm = E E E ^xHhj)pI^f,]p-: 

/=! i=l 3 = 1 

= y^ (f Vk'PmnxdT 


= Vk'PmrixdT, 


where we have used the fact that the boundary fluxes of adjacent elements cancel 
analytically. Thus, property 3 is satisfied. 

□ 

5. Results. The two-dimensional linear advection equation is used to verify and 
study the triangular-element SBP operators of Section 4. In particular, we consider 
the problem 


^ ^ ^ - n 

dt dx dy 
U{x, 0, t) = U{x, 1, t). 


U{x,y,d) 



V (x,y) efl = [0,1]^ 

and U(0,y,t) =U(l,y,t), 
1)5 ifr<i 
otherwise. 


where r(x, y) = y (x — + {y — The boundary conditions imply periodicity in 

both the X and y directions, and the PDF implies an advection velocity of (1,1). The 
initial condition is a continuous bump function that is periodic on El. 

A nonuniform mesh for the square domain El is generated, in order to eliminate 
possible error cancellations that may arise on uniform grids. The vertices of the mesh 
are given by 

x ^,3 = ^ {2ni/N) sin {271 j/N ), y,j = ^ sin {27ii/N) sin {27ij/N ), 

where N is the number of elements along an edge, and i,j = 0,1,2,N. The 
nominal element size is h = 1/A^. A triangular mesh is generated by dividing each 
quadrilateral {xij,Xi+ij,Xij+i,Xi+ij+i} along the diagonal from Xi+ij to Xij+i. 
Finally, for an SBP element of degree p, the reference-element nodes are mapped 
(affinely) to each triangle in the mesh. Figure 3 illustrates a representative mesh for 
p = 3 and N = 12. 

Remark 5. The use of an affine mapping on each element ensures that the 
mapping Jacobian is element-wise constant; consequently, the transformed PDE has 
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Table 3 

Maximally stable CFL numbers for the SBP operators on the nonuniform mesh with N = 32. 

p=l p—2 p=3 p—4 

C-SBP CFLmax 1-885 2.257 1.816 1.570 

D-SBP CFLmax 0.696 1.269 1.157 1.148 


constant coefficients in the reference space and the finite-difference operator remains 
SBP. More generally, curvilinear elements would require constructing SBP operators 
for each curved element. 

We consider both continuous (C-SBP) and discontinuous (D-SBP) discretizations 
using the SBP operators. The global operators for the C-SBP discretization are 
constructed using the assembly process described in Section 4.5. In addition, the 
periodic boundaries are transparent to the global C-SBP operator, that is, nodes that 
coincide on the periodic boundary are considered the same. The D-SBP discretization 
of the surface fluxes follows the nodal DC method outlined, for example, in Hesthaven 
and Warburton [13, Chapter 6]. We use SATs with a penalty parameter of cr = 1, 
which is equivalent to the use of upwind numerical flux functions. 

The classical 4th-order Runge-Kutta scheme is used to discretize the time deriva¬ 
tive. The maximally stable CFL number for each SBP operator was identified for 
N = 32 using Golden-Section optimization, where the CFL number was defined as 
•\/2Af/(/iAr) = At/Ar for a time-step size of At; Ar denotes the minimum dis¬ 

tance between nodes on a right triangle with vertices at (0, 0), (1,0) and (0,1). Each 
discretization was run for one period, t = T = 1 unit, and considered stable if the 
final norm of the solution was less than or equal to the initial solution norm. The 
results of the optimization are listed in Table 3. 

5.1. Accuracy and efficiency studies. Figures 5 and 6 plot the normalized 

error for the C-SBP and D-SBP discretizations, respectively, for SBP-operator 
degrees one to four and a range of N. Specifically, if Uq is the initial solution, u is 
the solution at t = T, and H is the appropriate SBP norm, then the error is 


Error = (u — H (m — Uq). 

The error is normalized by the integral norm of the initial condition. The mesh 
resolution ranges from iV = 4 to A = 64 in increments of 4. Each case was time 
marched using a CFL number of 0.9CFLmax, which was sufficiently small to produce 
negligible temporal discretization errors. 

The results in Figure 6 indicate that the D-SBP discretizations and the odd- 
order C-SBP discretizations exhibit asymptotic convergence rates of 0(hP~^^). In 
contrast, the C-SBP discretizations based on the p = 2 and p = 4 operators have 
convergence rates of 0{hP). These discretizations experience even-odd decoupling, 
i.e. checkerboarding, which is illustrated in Figure 4 by comparing the spatial error 
at t = T for the p = 1 and p = 2 C-SBP discretizations; the p = 1 error is smooth 
whereas the p = 2 error is oscillatory. Checkerboarding is not unique to the SBP 
simplex operators, and can be observed with one-dimensional SBP operators. We 
are currently investigating methods to address this issue with the even-order C-SBP 
discretizations. 

Figures 7 and 8 plot the normalized error versus CPU time for both the 
C-SBP and D-SBP discretizations. The runs were performed on an Intel® Core 
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Fig. 3. Example mesh with p = 3 and 
N = 12 for accuracy and energy-norm studies. 



Fig. 5. error between the C-SBP so¬ 
lution after one period, t = T, and the initial 
condition for different mesh spacing and oper¬ 
ators. 


Fig. 4. Solution error after one period, 
t = T, for the p = 1 (left) and p = 2 (right) 
C-SBP discretizations with = 12. 



Fig. 6. error between the D-SBP so¬ 
lution after one perior, t = T, and the initial 
condition for different mesh spacing and oper¬ 
ators. 


i5-3570K processor, and the code was implemented in Julia version 0.4.0. The 
efficiency provided by the high-order operators is apparent in both the continuous 
and discontinuous cases. 

5.2. Stability studies. Figure 9 shows the spectra of the C-SBP and SE spatial 
discretizations for the linear advection problem. Specifically, these eigenvalues are for 
the global operator Sx + Sy when N = 12. The eigenvalues of the C-SBP operators 
are imaginary to machine precision, which mimics the continuous spectrum for this 
periodic problem. This is as expected, because the boundary operators cancel 
between adjacent elements when the SBP derivative operator is assembled, leaving 
only the antisymmetric parts. The SE operator for p = 1 also has a purely imaginary 
spectrum, because it is identical to the linear SBP operator; however, the spectra of 
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Fig. 7. Normalized error of the C- 
SBP solutions after one period versus CPU 
time measured in seconds. 



CPU Time {s) 

Eig. 8. Normalized error of the D- 
SBP solutions after one period versus CPU 
time measured in seconds. 


the high-order SE operators have a real component. 

The consequences of the eigenvalue distributions are evident when the linear ad- 
vection problem is integrated for two periods. Figure 10 plots the difference between 
the solution norm at time t G [0,2] and the initial solution norm, i.e. the change 
in “energy”, 

AE = \-\Un — Mq Hmo, 

where u„ denotes the discrete solution at time step n. For this study, iV = 12 and 
the CFL number was fixed at 0.01 to reduce temporal errors. 

The energy history in Figure 10 clearly shows that the SE operators are unstable 
for this linear advection problem, while the SBP operators are stable. The small (lin¬ 
ear) decrease in the SBP energy error is due to temporal errors and can be eliminated 
by using a different time-marching method, e.g. leapfrog, or at the cost of using a 
sufficiently small CFL number. 

6. Conclusions. We proposed a definition for multi-dimensional SBP opera¬ 
tors that is a natural extension of one-dimensional SBP definitions. We studied the 
theoretical implications of the definition in the case of diagonal-norm operators, and 
showed that the multi-dimensional operators retain the attractive properties of tensor- 
product SBP operators. A significant theoretical result of this work is that, for a given 
domain, a cubature rule with positive weights and a full-rank Vandermonde matrix 
(evaluated at the cubature nodes) are necessary and sufficient for the existence of 
diagonal-norm SBP operators on that domain. We also developed simultaneous ap¬ 
proximation terms (SATs) for the weak imposition of boundary conditions, and we 
showed that an SBP-SAT discretization of the linear advection equation is time stable. 

We constructed diagonal-norm SBP operators for the triangle and tetrahedron. 
To the best of our knowledge, this is the first example of SBP operators of degree 
p > 2 on these domains. We also presented an assembly procedure that constructs 
SBP operators for a global domain from element-wise SBP operators. 

Finally, we verified the triangle-element SBP operators using both continuous 
and discontinuous (i.e. SAT) inter-element coupling. Results for linear advection on 
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Fig. 9. Eigenvalue distributions for the SBP (upper row) and the SE (lower row) spatial 
discretizations of the linear advection problem. Note the different ranges for the real axes. 



Fig. 10. Time history of the change in the solution energy norm. Note the use of a symmetric 
logarithmic scale on the vertical axis. 


a doubly periodic domain demonstrated the time stability and accuracy of the SBP 
discretizations. The results suggest that the proposed operators could be effective for 
the long-time simulation of turbulent flows on complex domains. 

Acknowledgments. All figures were produced using Matplotlib [16]. 

Appendix A. Decomposition of the SAT matrix /3a, + PyEy. 

Let E = /3a;Ea; -f- /3j,Ey. Recall the block-matrix dehnition of Ea, used in the proof 
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of Theorem 3.3. Using that definition, and a similar definition for E^, we have 


V-^EV-i =V-T 


E 

F G 


V-i = 


-T 


Ex FT 


V 


-1 


Byy 


-T 


E. FT- 
Fy Gy 


V' 


From the definition of E^, and Ey, we have that the entries in the block E are given by 


(e) = (f Vkix, y)Vm{x,y) {Bxrix + ByiT-y) dr V fc, m £ {1,..., 

\ /fe,m 

where Vk and Vi have total degrees less than or equal to r. We can decompose E by 
breaking the above integral into two integrals, one over r+ and one over r_: 


(e) = (/ 'Pk'Pm{Bxnx+Byny)dT + <L VkVmiBxnx + Byny) dT = E++ 

\ / k,m Jr_ 

y k,m e {1,..., n*}. 

where E+ and E_ are equated with the integrals over r+ and r_, respectively. 

Lemma A.l. The matrix E_ is negative semi-definite and the matrix E+ is 
positive semi-definite. 

Proof. We prove the result for E_, since the proof for the positive-definite matrix 
is analogous. Let u £ R"p be an arbitrary nonzero vector, and let Uk denote its 
entries. Then 

K K „ 

u^E-u = {ukVkix, y)) {uiViix, y)) {Bxnx + By^y) dF 

fe=l 1=1 •^T- 

= j) [U{x,y)f {Bxnx + Byny) dr 


where lA{x,y) = Y^^=i'^kPk{x,y). The integrand in the above is the product of 
a squared polynomial function and the non-positive quantity {Bxnx + Byny) < 0 
V (x, y) £ r_. Thus the desired result follows. □ 

We now turn to the main result of this appendix: 

Theorem A. 2. Suppose F^; = 0 and Ey = 0 in the definitions of Ex and Ey. 
Then, for any Bx,By S M, the matrix E = BxE.x + By^y can be decomposed into 
E = E+ -F E_ where E_|_ is positive semi-definite, E_ is negative semi-definite, and E± 
satisfy the accuracy conditions 

Pk^±Pm=(f VkVm {Bxnx P Byny) dr, V fc,m £ {1,2,... ,n*}. (A.l) 

Jt± 


Proof. Consider the decomposition 


E= E. 


E_ = V' 


E+ 

F-f 


(F+) 

G 


T-| 


V 


-1 


V' 


(F-) 

G-, 


T-| 




where E± are defined above, and the pairs (F+,F_) and (G-|_,G_) are yet-to-be- 
determined decompositions of BxEx + By^y and Bx^x + /3yGy, respectively; note that 
the above decomposition of E is distinct from its definition, which involves Ex and Ey. 
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The accuracy of the matrices E_|. and E_, as defined above, can be established 
using the same approach used in the proof of Theorem 3.3. Therefore, we focus on 
showing that E_ can be made negative semi-definite; again, an analogous proof can 
be used to show E_|_ is positive semi-definite. 

To prove that E_ is negative semi-definite, it suffices to show that the matrix 

E- (F_)^' 

F_ G_ 

is negative semi-definite. This will be the case if we can ensure G_ is negative-definite 
and the corresponding Schur complement, 

S_ = E_-(F_)^(G_)-'F_, 

is negative semi-definite; see, for example, [1, Appendix A.5.5]. 

We first tackle the definiteness of G„. Recall that G^, and Gy are symmetric but 
otherwise arbitrary. Therefore, the matrix G = P^G^ + PyGy has the eigendecomposi- 
tion G = RAR"'’, where R holds the eigenvectors and the diagonal matrix A holds the 
eigenvalues. For any set of eigenvalues, we can construct the nonunique decomposition 

G = RA+RT-1-RA_RT, 

such that A+ is diagonal positive-definite and A_ is diagonal negative-definite; note 
that any zero eigenvalue in A can be decomposed as c — c for arbitrary c > 0. Equat¬ 
ing G_ with RA_R"^ we have that G_ is symmetric negative-definite and therefore 
invertible. 

Finally, we need to show that S_ is negative semi-definite. From Lemma A.l, 
we have that E_ is negative semi-definite. Thus, showing S_ ^ 0 is equivalent to 
showing^ 

E_^(F_)T (G_)-'F_, 

This statement is true provided the entries in F_ are sufficiently small, which is 
certainly the case under the assumption that F^, = Fj^ = 0. This concludes the proof. 

□ 

Remark 6. In general, the assumption that F^, = Fy = 0 is stronger than 
necessary, since only S- ^0 is reguired; however, it is not clear how to weaken this 
assumption when j3x and jdy are not known a priori. 

Remark 7. The matrices E^, and Ey constructed for the simplex operators satisfy 
the conditions of the Theorem A. 2, i.e. Fa, = Fy = 0. 
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